Benjamin Planque - 2021
Institute of Marine Research. P.O, Box 6606 Langnes, 9296 Tromsoe, Norway
Tel: +47 48 89 30 43, email: benjamin.planque@hi.no, www.hi.no/en
RCaN is developped and maintained by Hilaire Drouineau and available from: https://github.com/inrae/RCaN
#install.packages("symengine")
#devtools::install_github("https://github.com/inrae/RCaN.git", subdir="RCaN")
#install.packages("ggplot2")
#install.packages("ggraph")
#install.packages("coda")
#install.packages("DT")
#install.packages("tidyverse")
require(symengine) # Symbolic math
## Loading required package: symengine
## SymEngine Version: 0.6.0
## _____ _____ _
## | __|_ _ _____| __|___ ___|_|___ ___
## |__ | | | | __| | . | | | -_|
## |_____|_ |_|_|_|_____|_|_|_ |_|_|_|___|
## |___| |___|
require(RCaN) # the RCaN library
## Loading required package: RCaN
## Loading required package: ROI
## ROI: R Optimization Infrastructure
## Registered solver plugins: nlminb, glpk, lpsolve.
## Default solver: auto.
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
require(ggplot2) # plots
## Loading required package: ggplot2
require(ggraph) # network plots
## Loading required package: ggraph
require(coda) # mcmc plotting
## Loading required package: coda
require(DT) # to visualise data tables
## Loading required package: DT
require(tidyverse) # get access to tidy syntax
## Loading required package: tidyverse
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble 3.1.4 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ✓ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::expand() masks symengine::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
This is a simple model with two populations. The first feeds on an external food source. The second population predates on the first one.
There are three explicit contraints:
1. the feeding rate should not exceed 100
2. the biomass of population1 should not exceed observations*1.2
3. the biomass of population1 should not be below observations/1.2
To begin with, only the first constraint is active
There several implicit constraints
1. the biomass of the populations should not go below the refuge biomass parameter
2. the feeding rates can not exceed satiation (sigma*biomass)
3. the relative change in population biomass from one year to the next is bounded by inertia (it cannot exceed exp(±mu))
rm(list=ls())
CaNfile <- "./TwoSpeciesRCaN.xlsx"
CaNmod <- buildCaN(CaNfile)
datatable(CaNmod$components_param[,1:8])
datatable(CaNmod$fluxes_def)
datatable(CaNmod$constraints)
datatable(CaNmod$series)
ggNetwork(CaNmod)
“polytope OK” indicates that this model has multiple solutions. This is the desired situation :-)
Other possibilities are
“empty polytope”
“polytope not bounded”
“unique solution”
“numerical error”
“potential problem”
checkPolytopeStatus(CaNmod)
## [1] "polytope ok"
gab <- getAllBoundsParam(CaNmod)
datatable(gab)
We use three sampling sequences (chains) and draw 1000 samples in each of them. The method used is the Gibbs sampler and there is no thining, i.e. all samples are kept.
CaNsample <- sampleCaN(CaNmod,1000, nchain = 3, ncore = 3, method = "gibbs")
Plot the time series for the population
Plot the time series for the feeding
Check the sampling performance: are the distributions convex? are the chains autocorrelated?
Comp2plot <- CaNmod$species
ggSeries(CaNsample,param = Comp2plot, plot_series = TRUE, ylab = 'Biomasses (1000t)')
Flux2plot <- CaNmod$fluxes_def$Flux
ggSeries(CaNsample,param = Flux2plot, plot_series = TRUE, ylab = 'fluxes (1000t/year)')
CanNames <- attributes(CaNsample$mcmc[[1]])$dimnames[2]
for(i in 1:4){
par(mfrow=c(2,2))
hist(CaNsample$mcmc[,i][[1]],breaks = 25,main = CanNames[[1]][i])
traceplot(CaNsample$mcmc[,i])
plot(autocorr(CaNsample$mcmc[,i:i][[1]],lags=1:25),type='l')
abline(h=0,col='red')
}
Everything looks great aside from the autocorrelation in the sampling. It might be worth “thining” the sampling, i.e. discarding samples.
Same sampling as above, but this time we apply a thining of 50, i.e. we only retain one every 50 samples.
CaNsampleThin <- sampleCaN(CaNmod,1000, nchain = 3, ncore = 3, method = "gibbs", thin = 50)
for(i in 1:4){
par(mfrow=c(2,2))
hist(CaNsampleThin$mcmc[,i][[1]],breaks = 25,main = CanNames[[1]][i])
traceplot(CaNsampleThin$mcmc[,i])
plot(autocorr(CaNsampleThin$mcmc[,i:i][[1]],lags=1:25),type='l')
abline(h=0,col='red')
}
We now include information for observational time series on the biomass of species 1. This is done by activating constraints C02 and C03 which state that the biomass trajectories sampled by RCaN should correspond to the observed biomasses /*1.2.
How does constraining the model on population 1 alters population 2 and the trophic flows?
CaNmod2 <- toggleConstraint(CaNmod,c("C02","C03"))
datatable(CaNmod2$constraints)
checkPolytopeStatus(CaNmod2)
gab <- getAllBoundsParam(CaNmod2)
datatable(gab[1:10,])
CaNsample2 <- sampleCaN(CaNmod2,1000, nchain = 3, ncore = 3, method = "gibbs", thin = 50)
ggSeries(CaNsample2,param = Comp2plot, plot_series = TRUE, ylab = 'Biomasses (1000t)')
ggSeries(CaNsample2,param = Flux2plot, plot_series = TRUE, ylab = 'fluxes (1000t/year)')